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Abstract 

We introduce a quasi-geostrophic model of core dynamics, which aims at describ- 
ing core processes on geomagnetic secular variation timescales. It extends the for- 
malism of Alfven torsional oscillations by incorporating non-zonal motions. Within 
this framework, the magnetohydrodynamics takes place in the equatorial plane; it 
involves quadratic magnetic quantities, which are averaged along the direction of ro- 
tation of the Earth. In addition, the equatorial flow is projected on the core-mantle 
boundary. It interacts with the magnetic field at the core surface, through the radial 
component of the magnetic induction equation. That part of the model connects 
the dynamics and the observed secular variation, with the radial component of the 
magnetic field acting as a passive tracer. We resort to variational data assimilation 
to construct formally the relationship between model predictions and observations. 
Variational data assimilation seeks to minimize an objective function, by computing 
its sensitivity to its control variables. The sensitivity is efficiently calculated after in- 
tegration of the adjoint model. We illustrate that framework with twin experiments, 
performed first in the case of the kinematic core flow inverse problem, and then in 
the case of Alfven torsional oscillations. In both cases, using the adjoint model allows 
us to retrieve core state variables which, while taking part in the dynamics, are not 
directly sampled at the core surface. We study the effect of several factors on the 
solution (width of the assimilation time window, amount and quality of data), and 
we discuss the potential of the model to deal with real geomagnetic observations. 

1 Introduction 

Current descriptions of core dynamics rely on two sources of information: observations 
of the magnetic field, and physical laws governing the evolution of the state of the core. 
The Earth's magnetic field is assumed to have an internal origin through the process of 



1 



geodynamo; it is generated and sustained by fluid motions in the metallic liquid outer core, 
and varies on a wide range of time scales reflecting the various time and space scales of 
core magnetohydrodynamics. 

The quality of observations of the Earth's magnetic field has much improved since the 
set-up of the first network of magnetic observatories by Gauss and co-workers in 1834, 
which was followed by the large increase in the number of observatories at the beginning 
of the twentieth century. Other turning points have occurred since: the introduction of 
the proton precession magnetometer, the development of declination/inclination magne- 
tometers (DIflux) widely used in observatories by the 1970's, and finally the rise of the 
Intermagnet network of digit al observatorie s shar ing modern measurement practices after 
1990 (see, e.g., the review by iTurner et al.l . 120071 ). The good temporal coverage of obser- 
vatory data has now been supplemented by the excellent spatial coverage of satellite data. 
Following the launch of three low Earth orbiting satellites -Oersted, CHAMP and SAC-C- 
supplying geomagnetic data, a continuous satellite time series extends now to 10 years. 

The magnetic field can be downward continued throughout the solid mantle to the fluid 
core surface. Field models are built, describing the ra dial component of the main field and 
its tim e evolution at the core-mantle boundary (CMB) f lHulot et aUl2007t I Jackson and Finlay . 



20071 ). Calculating the geomagnetic secular variation, the first time derivative of the main 



field time series, emphasizes rapid changes of the magnetic field, on characteristic time 
scales ranging from years to centuries. Inversions of a snapshot of the geomagnetic secular 
variation can be performed using the radial ind uction equation at the CMB, in order to re- 



trieve the large-scale part o f the flow beneath it ( lEvmin and Hulotl . 120051 ; iHolme and Olsen 
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rms) speed 
however, face non- 



of these flows is typic ally of the ord er on 15 km/y. Such inversions 
uniqueness problems f lBackud . Il968f ). Further assumptions are thus required to remove 
the non-uniqueness, and a great part o f the w ork consists in finding constraints and reg- 
ularizations to specify the flow (jHolmd . 120071 . section 8.04.2). Alongside these kinematic 
inversions, numerical model s of the geodynamo have been a vailable for more than 10 years, 
since the pioneering work of Idatzmaier and Roberts! ( 119951 ). The magnetic field generated 
by those dynamical models explains features of the Earth's magnetic field (dipolar ge- 
ometry, spatial spectrum); ye t their parameter s are far f r om th ose of the Earth ' s core 
riChristensen and Wichtl . 120071 . section 8.08.4). iRau et al.l fcOQOh and lAmit et al.l fl2007h 
tried to connect these two approaches (core flow inversion and forward numerical mod- 
elling) when inverting synthetic data from dynamo models. They found their core flow 
inversion method and the additional regularization to be adequate for the retrieval of 
large-scale flow and magnetic field patterns. 



A quality-control of core flow models is the angular momentum they carry (IJault et al 



19881 ). Comparison of these estimates with core angular momentum changes inferred from 
decadal length-of-day variations is encouraging, yet discrepancies remain. Angular mo- 
mentum series are also derived f rom atmospheric an d oceanic flow models, based upon the 
data assimilation methodology (IKalnav et al.l . Il996l ). Those variations of angular momen- 
tum account very well for the observed seasonal and interannual changes in length-of-day 
(ICherl l2005l : iGrossl . 120071 . section 3.09.4). 
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Data assimilation, routinely used in atmospheric science and more recently in oceanog- 
raphy, is now in early stages of use in the field of core physics. Applied to the core, this 
technique should allow us to interpret the secular variation in terms of dynamics, thereby 
enlarging the work d one on kinematic core flow inversion. Resorting to a toy model, 



Fournier et al.l ( 120071 ) assimilated synthetic data in a one-dimensional model that retains 



characteristic features of the induction and Navier-Stokes equations. They concluded that 
a good knowledge of the observed magnetic field can be translated into a good knowledge of 
core flow, through the process of data assimilation, which takes explicitly into account the 
dynamical relation ship that exists b etween magnetic and velocity fields. This conclusion 
was also drawn bv ISun et al.l ( 120071 ) . using a much similar toy model and a different im- 
pleme ntation of d ata as similation (sequential as opposed to variational). In a preliminary 
study, iLiu et al.l (120071 ) applied an optimal interpolation scheme to a three-dimensional 
dynamo model, using a synthetic set of observations. They showed in particular that 
assimilation of observations could partially alleviate the negative impact of wrong model 
parameter values. Still, the values of the parameters used typically in that class of simula- 
tions are far from being Earth-like, due to the numerical cost of their integration. There is 
hope, however, that systematic and appropriate exploration of the parameter space of those 



model s will eventually yield scaling laws of the kind proposed by IChristensen and Aubert 



( 120061 ). which will ultimately permit a reliable extrapolation between their output and the 
observed secular variation. In the context of geomagnetic data assimilation, the current 
situation is even worse, though, since one assimilation run requires several tens of forward 
realisations. 

A solution to this conundrum is to construct a simplified dynamical model, tailored to 
the study of the secular variation. In this paper, we introduce a model which relies on quasi- 
geostrophic dynamics. As a matter of fact, on rapid time scales as those characterizing 
the secular y ariation, rotational forces prevail on magnetic forces in the bulk of the fluid 
( IJaultl . 120081 ). The resulting two-dimensional flow interacts with the radial component of 
the magnetic field (in that instance, a passive tracer) at the core-mantle boundary. The 
quasi-geostroph ic assumption has recently been used in the framework of kinematic core 



flow inversions (IPais and Jaultl . l2008t iGillet et al.l . 120091 ). It also provides us with the tools 



necessary to build a dynamical model of the geomagnetic secular variation. 

Being quasi-geostrophic, this model comprises the equations for torsional Alfven waves, 
for which the dynamics is axi symmetric; Alfven torsional waves are associated with geostrophic 
(zonal) motions in the core (iBraginskyl . Il970l ). The frequency of these oscillations is pro- 
portional to the rms strength of the magnetic field B s perpendicular to the rotation axis 
s denotes the cylindr i cal radius). Acc ordingly, using a dat abase of computed core flows, 



Zatman and Bloxhaml (119971 . 119991 ) and iBuffett et al.l (120091 ) have calculated radial profiles 



of the quadratic cylindrical radial component of the magnetic field averaged on geostrophic 
cylinders, {B^ }, within the core. 

Our quasi-geostrophic model generalizes that axisymmetric approach by adding non- 
zonal motion and magnetic field. Theoretical solutions of the model include various fami lies 
of diffusionless hydromagnetic waves, some of which were first studied by iHidd ( 119661 ) in 
order to explain the observed secular variation. 
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Figure 1: Geometry of the system and notations; a: Side view, b: Equatorial section. £ is 
the equatorial plane and the CMB is the outer sphere, located at r = r Q . is the radius 
of the inner core. Spherical (O, e r , eg, e^,) and cylindrical (O, e s ,e v ,e z ) coordinate system 
bases are defined at the core surface and in S respectively. 

The goal of this paper is thus to describe a quasi-geostrophic forward model of the 
Earth's core fast dynamics, and to place it at the heart of a geomagnetic data assimilation 
process. In section 2, we derive that quasi-geostrophic model, along with its link to the 
observations at the CMB. Variational data assimilation is introduced in section 3 and its 
principles are illustrated in section 4 with twin experiments. That section begins with the 
study of the classical kinematic inversion of a steady core flow, set within that framework. 
A second illustration is dedicated to the retrieval of the magnetic field sheared by Alfven 
torsional waves. Results are summarized and discussed in section 5. 



2 Quasi-geostrophic forward model 

We shall model the Earth's outer core as a spherical fluid shell of inner radius r, and 
outer radius r . The fluid has density p, and it is electrically conducting. Its magnetic 
diffusivity is rj. The system is rotating at angular velocity Q around the z-axis. Figure [T] 
sketches the geometry and summarizes the notations. 

The fluid is characterized by its magnetic field B, its velocity u, and the reduced 
pressure n, that includes pressure and the centrifugal potential. We choose r Q as length 
scale and Bq, a typical magnetic field intensity in the core interior, as the magnetic field 
scale. Velocities are scaled with the Alfven waves speed 

v A a 



in which p is the magnetic permeability of free space. The pressure scale is pV%, and the 
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time scale is the Alfven waves period: Ta = t /Va- 

The evolution of the magnetic field in the core is governed by the induction equation, 
which has the dimensionless form, 



d t B = Vx (u x B) + S~ L VB. 



(2) 



The Lundquist number S charac terizes the rati o between the magnetic diffusion time 
and the period of Alfven waves (e.g. iRobertsl . 119671 ). 



S 



(3) 



For the Earth's core, S is of the order of 10 4 to 5 x 10 4 (IJauftl . 120081 ). On secular variation 
time scales, diffusion becomes negligible compared to induction, hence the large value of 
S, which yields the frozen-flux approximation. 

Assuming that the mantle is electrically insulating on secular variation time scales, the 
magnetic field can be downward continued to the core-mantle boundary. In the case of a 
perfectly conducting fluid (the frozen-flux limit), the radial component of the magnetic field 
B r is the sole magnetic component continuous across the spherical CMB. At the top of the 
core, B r interacts with core motions by means of the radial component of the diffusionless 
version of equation (J2]) at r = r D , 



d t B r 



H 



[\lB r 



with the horizontal divergence operator Vjt defined as 

Vh ■ v = (sin 6y l d e (sin 6v e ) + (sin 6)' 1 d^v^, 



(4) 



(5) 



where (r, 9, ip) are the spherical coordinates. It is that equation at the core-mantle bound- 
ary that connects our model to the observations. The time-varying B r acts as a passive 
tracer (a drifting buoy), because it interacts with the velocity field at the core surface and 
does not affect the dynamics that sets up in the interior of the core (see below). 

On secular variation time scales rotation forces are much larger than magnetic forces 
in the bulk of the fluid. IJaultl ( 120081 ) suggests that rapidly rotating motions of lengthscale 
I are axially invariant if the non-dimensional Lehnert number, A;, is small enough. That 
number measures t he ratio betwee n the period of inertial waves, 1/fi, and the period of 
Alfven waves, I/Va fjLehnertl . 11954 ): 



A, 



Bn 



fi(/i p) 1/2 / 



A 



(6) 



Note that A; is a decreasing function of /. In his calculations, the flow appears to be 
invariant in the direction parallel to the rotation axis, provid ed A; 1 . For the Earth's 



core, with Bq of the order of 2 mT fjChristensen et all 120091 ) and I ~ 10 6 m, A 



10' 



Therefore, we shall assume that the flow is geostrophic at leading order. Working in the 
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equatorial plane E (crosshatched in Figured]), a cylindrical set of coordinates (s, <p, z), with 
e z parallel to the axis of rotation, is well-suited to study the resulting columnar patterns. 
The main force balance involves the Coriolis force and the pressure gradient 

2e z x u° = -Vn°, (7) 

where the superscript denotes the main order. Taking the curl of equation ([7]) yields the 
Proudman- Taylor theorem, namely the z-invariance of the flow. 

Within a spherical container, u° does not satisfy the non-penetration boundary condi- 
tion at the CMB, except if it consists of cylindrical flows organized around the rotation 
axis. Thus we have to add the first-order contribution in A ro of the Coriolis force, leading 
to 

D t u° + 2A r T o 1 e, x u 1 = - VII 1 + (V x B) x B, (8) 

where D t denotes the material derivative D t = d t + (u° • V). At first-order, magnetic forces 
are a natural candidate to trigger a departure from geostrophy, since magnetic energy is 
large compared to kinetic energy in Earth's core. Buoyancy forces are another candidate 
that we could additionally take into account, which we discard for now for the sake of 
simplicity. Viscous forces are neglected, while equation OH]) shows that the Coriolis force 
is scaled with the inverse of the Lehnert number A ro . The non-penetration boundary 
condition at the CMB: u 1 • e r = at z = ±h, yields a linear dependence of u 1 with respect 
to z 

u\(s > r i: if, z) = z/3u° s (s, tp). (9) 

If h = yrf — s 2 denotes the half-height of the column located at a given cylindrical radius 
s, the slope of the upper surface is dh/ds, and we can define 

P(s) = h~ l dh/ds. (10) 

The notation /3 has been chosen in reference to the /3-plane approximation. This approx- 



imation -with uniform (3- is ubiquitous in geophysical fluid dynamics (e.g. iValfisl . 12006 . 
section 2.3). It is convenient, indeed, to study planetary Rossby waves assuming that the 
Coriolis parameter (/o = 2f2cos#) varies linearly with latitude; (3 is then the northward 
gradient of the Coriolis parameter. 

According to our quasi-geostrophic approach, the flow in the outer core is nearly two- 
dimensional, which makes it natural to take the vertical average of the Navier-Stokes 
equation flS]). The vertical average (•) of a quantity X is defined as 

(X) (s, cp) = -i- / X (s, tp, z) dz. (11) 



2h(s) 



-h 



In a multiply-connected domain, the (^-averaged y orticity equation is not equivalent to 



the (/9-averaged Navier-Stokes equation (IPlautl . 120031 ) . as the former does not ensure the 
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existence of the pressure field; accordingly, we describe the evolution of the non- zonal flow 
u NZ by means of the axial vorticity equation, while the yj-averaged momentum equation 
directly provides us with the time changes of the zonal velocity u z = u z e v . In the re- 
mainder of this paper, the superscript capital Z marks zonal quantities. It should not be 
confused with small z, which refers to the direction of rotation. 

The non-zonal (NZ) velocity field u NZ is written as the curl of a z-independent non- 
zonal streamfunction ty, 

u NZ (s, cp) = V x <p)e z . (12) 

The non- zonal vorticity field £ is defined by ( = V x u NZ , and its vertical component is 

( Z (s,<p) = -Vjs*{s,<p), (13) 

in which the equatorial Laplacian operator is defined by 

V| = (sd s ) + s~X- (14) 

If we now curl the non- zonal part of the z-averaged Navier-Stokes equation (JHj), we find 
that the vertic a l com ponent of the vorticity equation is then identical to equation (17) of 



Pais and Jaultl (120081 ). with an explicit right-hand side term, 

D t ( z - 2A r ->- 1 ^\f = (s^d s d v + s~%) ({Bl) - (B 2 )) 

+ (?,s- 1 d s - s~ 2 dl + d 2 ) (B S B V ) . (15) 

The magnetic surface terms, which appear when taking the z-average of the Lorentz 
force, are neglected because we assume the magnetic field at the core surface to be much 
smaller than in the bulk of the fluid. The non-penetration boundary condition, at s = r Q 
and at the tangent cylinder s = r iy imposes \& = at both boundaries. 

The time evolution of the zonal velocity u z (s) = su g (s) is governed by 

D t u g = (s'h)' 1 d s (s 2 h (B S B V )) . (16) 

The two flow equations (fT5l) and ([TBI) contain z-averaged squared magnetic quantities 
(B 2 ), (B 2 ~) and {BgB^), whose time evolution is in turn derived from the diffusionless 
version of the induction equation ([2]) 

d t (B 2 ) = -[u°.V](B 2 ) + 2(B 2 )d s u s + 2s-yB s B ip )dy s} (17) 
d t {Bl) = -[^■V]{Bl)-2{Bl)d s u» + 2s(B s B^)d s (s- l ul), (18) 
dt(B a B v ) = ~[u -V}(B s B v ) + s(B 2 )d s (s- 1 ul)+s- 1 (B 2 )d v u° s , (19) 

where we have made use of the solenoidal character of B and u. 

The vertical averaging naturally sets the magnetohydrodynamics in the equatorial plane 
£ (Figure [Tb). The flow is then projected at the CMB, where it interacts with the radial 
magnetic field B r via equation (j3J) above. 

An alternative model, where the velocity field entering the set of equations (ITT]) to (TT9"]) 
has a ^-component given by ([9]) and a ^-component modified in order to ensure that the 
total velocity field remains solenoidal, is discussed in appendix [S] 
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3 Variational data assimilation framework 



In this section, we introduce the geomagnetic sec ular variation data assimilation prob- 



F 

lem with the notations suggested by llde et al.l (119971). In comparison with the 4DVar label 
commonly used in data assimilation (e.g. ICourtierl . 119971 ). our framework is l+2DVar, since 
state variables are defined in two-dimensional spaces to which a third (temporal) dimension 
is added - the 3DVar labe l is traditionally reserved for three-dimensional (in space) static 
problems flCourtierl . [l997h . 

The state vector x for the Earth's core gathers the variables involved in the description 
of the system state 



(20) 



where superscript T means transpose. Observations y are available at T y different epochs 
and N y spatial locations during the assimilation time window [0, T]; the size of y is typically 
smaller than the dimension of x. The observation vector is related to the true core state 
x* via the observation operator H: 



Hx* + e, 



(21) 



in which e is the observation error. 

Variational da ta assimilation aims at adjusting a model solution x-^ to the observations 
( iTalagrandl . 119971 ). by minimizing a misfit function which comprises the quadratic dis- 
crepancy -if H is linear and errors are Gaussian- between the true observations and those 



predicted by the computed state, Jh (jCourtierl . 119971 ): 



J, 



H 



3=1 



H x J 

n 3*j 



(22) 



where j is the discrete time index and R = E (ee T ) is the observation error covariance 
matrix, E(-) denoting statistical expectation. The matrix R describes the level of confidence 
we set in the observations. 

It might be necessary to constrain the solution sought in order to enforce its uniqueness, 
especially if the problem is non-linear. Constraining the assimilation refers to either adding 
a background state x 6 from which the estimate shall not strongly deviate, or applying 
additional constraints on the core state. Imposin g a penalty, the goal of which i s to favor a 
solution with a moderate level of complexity (e.g. lCourtier and Talagrandl . ll987l ). described 
by a matrix C applied to the state vector, consists in adding a second term to the objective 
function, of the form 



J, 



c 



3=0 



(23) 
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where W = C T C. The total misfit function J to minimize is then given by the sum 



J ~ 2 H + 2 C ' 



(24) 



in wh ich the two contributions are weighted by two scalars, an and ac ( IFournier et all 
2007h . 

In practice, x-^ is the solution of a numerical non-linear model M, that describes the 
temporal evolution of x at any discrete time tj G [0,T], 



x J 



M,x^; 



(25) 



that notation symbolically summarizes the set of equations developed in section [2j Since the 
temporal history of the state is constructed with a dynamical model that couples its various 
components, we can, in principle, retrieve information about every state variable, even if 



only part of the state vector is directly probed by the observations on hand ( IFournier et al. 



20071 ). Moreover, the initial condition xo, termed the control vector, uniquely defines the 
model trajectory in state space. That reduces the dimension of the associated inverse 
problem: we only seek the initial condition, xo, starting from which the temporal evolu- 
tion x a will best fit the observations; in assimilation parlance, this best solution is called 
the analysis. The minimization of J (that is the search for x°) is performed with a descent 
algorithm that involves the sensitivity of J to its control variable s, x n : V yn J. Its transpose 



is effi ciently estimated with the use of the adjoint model M T ( )Le Dimet and Talagrandl . 
19861 ). For a given xo, one couples a forward integration of M with a backward integra- 
tion of M T to express t he gradient of J. The adioint model is that of the lo c al tan gent 
linear equations (see e.g. iTalagrand and Courtierl . 119871 ; iGiering and Kaminskil . Il998l ). In- 
troducing the adjo int variable a of x, t he adjoint equation imposed by the cost function 



(equation (pift ) is (IFournier et all 120071 ) 



a/_i = MjLxBj + a H H]L J 'R j lx (H^x^j - y^) + acW^jX^j. 



(26) 

where H T is the transpose of the observation operator (equation (I2ip ). which projects a 
vector from observation space to state space. Through equation (|26i) . the adjoint field is 
fed by observation residuals (the difference Hx — y), as soon as there is an observation 
available. 

The backward integration requires the knowledge of the model trajectory over the 
assimilation time window. The storage of the complete trajectory may cause memory 
issues, which are traditionally resolved using a checkpointing strategy. The state of the 
system is stored at a limited number of discrete times, termed checkpoints. Over the 
course of the backward integration of the adjoint model, these checkpoints are used to 
reco mpute local port ions of the trajectory on-the-fly, whenever those portions are needed 



e.g. iHersbachl . |l998). 



Non-linear optimization is performed with the M1QN3 routine ([Gilbert and Lemarechal 



19891 ). which implements a limited memory quasi- Newton algorithm that approximates the 
inverse Hessian (second derivative) of J. 
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4 Applications 



We show two illustrations of variational data assimilation applied to fast core dynamics, 
as described by our quasi-geostrophic model, with synthetic data. The methodology of twin 
experiments is explained and applied to a steady non-zonal flow problem and, in a second 
step, to a dynamical zonal model of torsional Alfven waves. These two problems correspond 
to two subsets of the model presented in section [2J 



4.1 Twin experiments 

Instead of being satellite or observatory data, observations in our twin experiments are 
created from a synthetic true state, which is the result of the integration of the forward 
model for a given set of initial conditions, x . Synthetic data have the advantage of 
representing only the physics involved in the model and are, in a first step, appropriate 
to test the implementation of the variational data assimilation algorithm. A database of 
observations is produced with equation (|2"T|) . 

To construct the observation catalog, we include some geophysical realism by averag- 
ing the state at the core-mantle boundary. We choose the averaging window so that it 
corresponds to the ignorance of the spherical harmonic coefficients of degree n > L. Then, 
the product Hx* corresponds to the convolut ion over the core-ma ntle boundary of the true 



state with a Jacobi polynomial of degree L (IBackus et al.1 . Il996l . paragraph 4.4.4): 



(Hx*)(0 o ,p o ) = / / x*(0,p)Pi 1,o) (cosa)sin0d0(fy>, (27) 

4vr J e=0 J^q 

in which (0 O , y? ) are the coordinates of the observation locations, a is the angular distance 



betw een the points (0 O , ip Q ) and (0, (p), and p£ is a Jacobi polynomial (lAbramowitz and Stegun 



1964L p. 773). In the following experiments, we set L = 15 and observations are made at 
a fixed temporal frequency. 

Next, we start the assimilation with a different set of initial conditions, called initial 
guess, Xq. After a forward integration, the computed observable, HB r in section I4T21 and 
H.B/ in section 14. 3[ is compared (over the entire time window) with the observations, 
and the discrepancy between the two gives the initial misfit (see equation (124")) ). After 
assimilation, the decrease of the misfit and the relative difference between x* and x a , in 
the / 2 -sense, are used to assess the quality of the recovered state. 



4.2 The kinematic core flow problem 

In this section, a connection is made between core flow kinemati c inversions and data 
assimi latio n. The steady flow hypot hesis has been previously used by IVoorhies and Backus 
( 119851 ) and lWaddington et al.l ( 119951 ) to break the non- uniqueness of the kinematic inversion 
problem. Here, we study the effect of a steady non-zonal and equatorially symmetric flow 
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on the evolution of the radial magnetic field, and more particularly its secular variation 



B r , 



B r = -V H • (u NZ B r ) . (2* 



Symmetry with respect to the equator ensures uniqueness of the solution when B r and B r 
are perfectly known. The time scale characterizing this problem is the advection time, t a d v - 
The non-zonal velocity effectively enters equation (128]) via the non-zonal streamfunction 
*lf(s,(p) (see equation (fl2l ). projected at the core-mantle boundary. The state vector x 
includes a parameter, the streamfunction and a variable, the radial magnetic field 
B r . Here, \I/ is called a parameter because it is steady in that experiment. We seek 
the distribution of \I/(s, ip) which best explains the observed synthetic database of secular 
variation B° at the top of the core. 
The tangent linear form of (1281) is 



5B r = P (6V,B r ) + Q (V,6B r ) , (29) 
d t 5B r = 5B r , (30) 

where 5^, 5B r and 5B r are the differentials of B r and B r , respectively, and P and Q 
the differentials of the right-hand side term of equation (|28[) with respect to \I/ and B r , 

respectively (they are developed in appendix [U]) . Let us introduce fy T ,Bj and B r as the 
adjoint variables of B r and B r . The adjoint model of equations (T2"§j) and (1311 is then 



E pT (^ T >^)> ( 31 ) 

3=0 

B T r = Q T (B r T ,v), (32) 

in which P T , Q T are the adjoint functions of P and Q (see detailed equations in appendixO). 
The link to the observations has been obtained from equation (1241) . and is computed at 
each time step as in equation (1261) : 



B r = a H (H 1 HB r - H 1 B°\ . (33) 
The trajectory of the true state is computed from the following set of initial conditions: 



1. B*(6,ip,t = 0) is obtained from the CHAOS main field model (lOlsen et all 120061 ) 
for epoch 2002, truncated at spherical harmonic degree and order 12. It is taken 
outside the tangent cylinder and multiplied by a sine function of 9 in order to have 
Bl(9,ip,t = 0) = at the tangent cylinder (see Figure EJ, 

2. \l/*(s,(/?) is show n in Figure [3k, it is the non-zonal part of an inverted flow from 
Pais et al.l (12004 ) truncated at degree and order 4, and multiplied by cos 2 9 = (1 — s 2 ) 



and a function of s, (s — r,), in order to satisfy the flow boundary conditions at 
s = Ti, r a . It is norrr. 
the scaling V a dv is su 
with 9 C = asin (rj/r Q ) 



s = Ti,r . It is normalized in order to have a dimensionless rms velocity of order 1 
the scaling V a dv is such that J Q n c {u 2 s + u 2 ^) sin 9d9dip = V 2 dv f Q n f£ c sin 9d9df 
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Figure 2: Map of the main field at epoch 20 02 truncated at sph erical harmonic degree and 
order 12 modified from the CHAOS model flOlsen et all 120061 ) . True radial component of 
the magnetic field B\. at initial time for the twin experiments. Contours are drawn each 
0.5B r r ms (solid (resp. dashed) for positive (resp. negative) values). Note that the problem 
is linear in B r . 



We consider perfect observations, setting e = in equation ( |2TI) . For the following 
simulations, the numerical time step is 6 x 10~ 5 t a d v for integration times ranging from 
0.03 t a d v to 0.57 t a d v . Other numerical parameters relative to the simulations are given in 
appendix [B] 

In this problem of seeking a steady streamfunction that explains the observations, we 
want to show the benefit of including the temporal dimension (data assimilation) instead 
of relying on a single observation epoch, as it is the case for a standard kinematic inversion. 

We first consider solutions obtained with only one observation epoch and less obser- 
vation locations (N$ = 50, N° = 11) than grid points (N g = 200, iV^ = 15). We start 
the assimilation with a first guess \1/ 9 corresponding to the minimal hypothesis: ^ g = 0. 
The solution we obtain gives us some information about what could be achieved within 
the kinematic framework in that configuration. Figure [3b shows in particular that the true 
state (Figure Ek) is not completely recovered, due to the truncation used in the construc- 
tion of H and the limited number of observation locations, compared to the total number 
of grid points. 

Using this solution obtained from a single epoch inversion as a reference solution, we 
can study the benefit we get resorting to a timeseries of observations, as opposed to a single 
snapshot. In other words, we investigate whether the issue of spatial subsampling can be 
partially fixed by considering the temporal dimension. To that end, we do experiments 
with assimilation time windows ranging from 0.03 t a d v to 0.57 t a d v , at a given temporal 
frequency of observation f y = 100 t~^ v , keeping the same number of virtual magnetometers. 

Results of a typical experiment are shown in Figure |3b, for which T = 0.57 t a d v - The 
large scale pattern is retrieved, but the solution is polluted by small spatial scales (no extra 
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Figure 3: Maps in the equatorial plane £ of the steady streamfunction a: true state, 
b: analyzed state with a single epoch inversion, c: analyzed state with T = 0.57 t a d v - 
Contours are drawn each 0.05 (solid (resp. dashed) for positive (resp. negative) values). 
Extrema are -0.33; 0.54 (a), -0.49; 0.54 (b) and -0.43; 0.56 (c). 



smoothing term is added to the misfit function). We find, however, that the consideration 
of the temporal dimension improves the solution. Moreover, the distance between the true 
streamfunction and the retrieved streamfunction becomes smaller when the assimilation 
time window is wid ened (see Figur e HI) . 

As described by lEvensen (2007. chapter 6), if one enlarges the width of the assimilation 
time window in a non-linear context, the misfit function presents more and more spikes 
and minima. A very good first guess is thus needed to converge to the global minimum. 
To circumvent this issue, we decided to use the results obtained over short assimilation 
time windows as initial guesses for assimilation over longer time windows. That strategy 
is analogous to the approach used in the atmospheric variational assimilation community, 
which consists in solving a series o f stron g constraint inverse problems, defined for separate 
subintervals in time (e.g. lEvensen! . 120071 ). 



4.3 Forward and adjoint modeling of Alfven torsional waves 

In non-rotating magnetized flows, classical Alfven waves result from the balance be- 
tween inertial and magnetic forces. 



In the Earth's core, where the Coriolis force is large, iBraginskyl (Il970l ) showed that a 
special class of Alfven waves comes into play, in which only the component of the magnetic 
field normal to the axis of rotation, B s , participates. Associated motions are geostrophic; 
they are organized in axial cylinders about the axis of rotation, hence the name torsional 
oscillations. The period of torsional waves depends on the strength and distribution of B s 
inside the core. 

In order to study these waves, one can consider a subset of the complete dynamical 
model. Since torsional waves are geostrophic and axisymmetric motions, let us discard the 



non- zonal part of the flow and magnetic induction in equations (jTojl to ([T9 
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Figure 4: Effect of the assimilation time T on the analyzed state at a fixed observation 
frequency f y = 100 t~^ v , measured with the norm of the relative error ||^' — *& a \\ 2 /||^E r *|| 2 . 
This norm is equal to one when the initial guess ^ 9 = 0. The simulation at T = is 
referred in the text as the single epoch inversion. 

In addition to the vertical average, (•), introduced in equation (ITT]) , we now define the 
average on a geostrophic cylinder, {•}, by 



2tt 



2tt 



(X)(8,<p)d<p. 



(34) 



In the following, we will not indicate explicitly the dependence on s of quantities averaged 
on a geostrophic cylinder. Application of {•} to equation ( fi~6j) yields 

d t oo g = (s 3 h) - 1 d s (s 2 h {B 8 B V }) . (35) 
Similarly, the equations governing the evolution of magnetic quantities become 



dt{Bl) = 0, 
d t {B s B v } = s{Bl}d t 



UJg. 



(36) 
(37) 



Written in terms of the geostrophic angular velocity u g , the torsional wave equation is 

d 2 t u g = (s 3 h)- 1 d s (s 3 h{B 2 s }d s u g ). (3* 
Equation fl38|) can be transformed into a set of two first-order equations 



d t ujg = (s 3 h) d s r, 
d t r = s 3 h{B 2 s }d s u 



(39) 
(40) 



in which r = s 2 h {BsB^,} is an auxiliary variable. 
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We have taken as boundary condition for the angular velocity: 



d s u g = 0, at s = r a . (41) 

Then, the boundary condition, 

d s uj g = 0, at s = Tj, (42) 

ensures the conservation of the total angular momentum carried by the fluid in the com- 
putational domain. 

Projected at the CMB, those motions interact with B r through equation (j2J), which 
simplifies here into 

d t B r = -u g d v B r . (43) 

The system state gathers the geostrophic angular velocity oj g , the variable r, the cylindrical 
average of the s-component of the magnetic field, {B 2 }, and the radial component of the 
magnetic field B r at the CMB. 

We define ', r T , {B 2 } T , Bj as the adjoint variables of cu g , r, {B 2 S } , B r respectively. 
The torsional oscillations adjoint model is 

-ftof = dJ[s"h{B 2 }r T ]-Bjd v B r , (44) 
-d t r T = dj^h)- 1 ^], (45) 

T 



{B 2 S } (s) = ^fhid^if + a c W{B*}, (46) 

3=0 

-d t Bj = -d^[u 9 Bj], (47) 

where dj and are the adjoints of the differential operators d s and d v (see appendix [D] 
for more details on the adjoint system). The link to the constraint on {B^} has been 
obtained from equation (124]) . The model is completed by the information supplied by 
the observations, as in equation f)33jl . The boundary conditions for the adjoint model are 
t t = 0, at both s = Ti and s = r D . The temporal and spatial discretizations of this problem 
are described in appendix IB1 

For the experiments that follow, the set of initial profiles, which define the true state, 

is: 

1. the same Bl(9, <p,t = 0) as in the kinematic core flow problem of section I4T21 

2. a Gaussian function for the angular velocities: u;*(s, 0) = uo exp [— a~ 2 (s — s w ) 2 ], 
with cr~ 2 = 150 and s w = 0.65; it satisfies the boundary conditions fj4T]) and fj42|) . its 
amplitude being scaled by u (discussed hereafter), 

3. r*(s,0) = 0, 
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Figure 5: Torsional wave twin experiments results: successive profiles of the angular 
velocity u g , showing the propagation of a torsional wave in the computational domain 
s G [0.35, 1] during 1.16 Alfven time T4. Initial condition (black curve) for u g , and snap- 
shots, at 0.12 (green), 0.18 (blue) and 1.16 T4 (orange). 



4. an arbitrary function {B^Y (see the black curve in Figure El right) given by 

{B^Y (s) = ci + c 2 sin (vr/2 - L) + c 3 exp [-a B 2 (s - s B f] , with c x = 0.1, c 2 = 0.02 



C3 = 1, a B = 20, sb = 0.8 and L = 14s. It is normalized in order to have a 
dimensionless rms magnetic field of unity inside the core; the scaling B is such that 
C C B* a dsd<p = Bl C j;° dsd^. 

As the velocity has been scaled by the Alfven velocity, Uq is the ratio between the 
Alfven time and the advection time, u = TA/t a d v - 

Assimilation is performed on both u g and {B 2 }. We seek the steady profile of {B 2 } and 
the initial profile of u g which best explain the synthetic database of B°. Our first guess 
consists of flat profiles: u^(s,0) = O.luo and {B 2 } 9 = 0.6 (see the red curves in Figure[6]). 

We show here experiments with a fixed frequency of observations f y = 20 and as 
many observations locations as grid points. The observations are blurred by the averaging 
kernel (equation (I2T]) ). which causes errors. Consequently, the analysis can develop small 
scales, which are not very well constrained by the observations. We choose to reduce 
the complexity of the solution by adding a smoothing term to the cost function, taking 
etc = 10~ 8 in equation fl24l) . Here we penalize only the strong spatial gradients of {-£>;?}• 
The reference case (with uq = 0.34 and T = 1.16 T4, see Figures [6] and [7]) shows that both 
the angular velocity and interior magnetic field are well recovered. As shown in Figure 
the error field is substantially weaker after assimi lation. 

On a technical note, the M1QN3 algorithm ( Gilbert and Lemarechall . 1989 ). used in 



the optimization loop, stops in that case when the initial misfit is divided by a factor of 
4 x 10 5 , which is reached in 214 iterations. 

In order to assess the effect of the width of the assimilation window on the retrieved 
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Figure 6: Torsional wave twin experiments results: true state (black), profile before as- 
similation (red) and solution after assimilation (dashed-orange) for u g (left) and {B^} (s) 
(right). The parameters for that reference case are Uq = 0.34 and T = 1.16 Regular- 
ization has been added to the spatial derivative of {B^}, of amplitude ac = 10~ 8 (See text 
for details). 




Figure 7: Relative difference between observed and computed HB r at final time, 
[B°(T) — H£>/(T)] / ||£>°(T)|| 2 before assimilation (a) and after assimilation (b) for the 
reference case (same parameters as in Figure E]). Contours are drawn each 0.1 (a) and 10 _1 
(b) (solid (resp. dashed) for positive (resp. negative) values). Extrema are —0.94; 0.86 (a) 
and [-6.4;7.2]xl0~ 4 (b). 
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Figure 8: Effect on the analysis of the assimilation time window of width T (left) and 
of the amplitude of the initial pulse uo (right). Left: true state (black), profile before 
assimilation (red) and solution after assimilation (dashed-green, dashed-double-dotted- 
blue, dashed-dotted-orange) for {B 2 } (s), for different values of T (0.12,0.18 and 1.16 Ta 
respectively) for a fixed u = 0.34. Right: true state (black), profile before assimilation 
(red) and solution after assimilation (violet) for {B 2 S } (s). Dotted, dashed and dashed- 
dotted violet curves are obtained for different values of the amplitude of the initial pulse 
uo (0.014,0.14 and 1.4 respectively), keeping T fixed, equal to 0.18 Ta- 



state variables, we vary the assimilation time T between 0.12 and 1.16 T4, keeping w 
constant, equal to 0.34 as above. The geostrophic angular velocity is in all cases completely 
recovered (not shown) with similar spurious oscillations as in Figure |6] (left) near the outer 
boundary. On the other hand, the area over which {B 2 } is correctly retrieved increases 
with T, indicating that the assimilated area is controlled by the distance over which the 
initial pulse has propagated (Figure [HI left). Figure shows that, if T = 1.16 Ta, the 
wave has enough time to explore the whole domain. On the contrary, if T = 0.12 or 
0.18 Ta, a lesser portion of the domain is sampled by the wave, over which {B 2 } has 
been effectively retrieved. The angular velocity is better recovered than {B 2 } because 
it is directly connected to the observations, as opposed to {B 2 }. That can be seen in 
the adjoint equations: uj (equation f|44|) ) depends on Bj that contains B°, the observed 

quantity, whereas {B 2 } T (equation (}4"6"]) ) is only directly connected to r T . In turn, r T sees 
cuj, which is ultimately linked to the observations. 

For a fixed T, the dependence on the amplitude u has been studied (Figure El right). 
For uo = 0.014, u g and {B 2 } are not well recovered. Starting from uq = 0.14, however, 
increasing ujq by more than one order of magnitude has little effect on the retrieval of {B 2 } 
and no effect at all on bj a g . 

Let us stress that the convergence of the calculations presented here is also sensitive to 
the profile of the initial condition (for both true state and first guess), and to the amount 
of measurements, as in the steady case of section 14.21 Moreover, we have observed in other 
instances (not shown) that convergence is sped up if the acceleration r is not zero at initial 
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Figure 9: Torsional wave twin experiments results obtained with noisy data: True profiles 
(black), profile before assimilation (red) and solution after assimilation (orange, green) for 
uj g (left) and {B 2 } (right). Dashed-orange curves have been computed with perfect obser- 
vations, dashed-dotted-green curves considering centered, normally distributed Gaussian 
observation errors of standard deviation 10^ 1 B^. ms added to the database. Regularization 
has been added on the spatial derivative of }, of amplitude «c = 10~ 8 . The assimilation 
window width is T = 1.16 Ta, and uq = 0.34. 



time. More generally, that particular example shows that the success of assimilation is 
controlled by the intrinsic dynamics of the system under study, as well as a good guess of 
its state. 

Until now, we have assumed perfect observations: e = in equation (l2Tj) . In the 
prospect of future applications, observation error should be considered. In the next ex- 
periments, observations contaminated by errors are assimilated. Centered, normally dis- 
tributed Gaussian observation errors of standard deviation 10~ 1 Bj, ms are added to the 
previous database. That experiment is carried out with an assimilation window width of 
1.16 Ta- Even with a database contaminated by observation errors, it is still possible to 
recover the shape and strength of the true state (see Figure {B 2 S } seems less sensitive 
to observation errors than u g : we still have an extra penalty term in the misfit function as 
in the perfect observations case, but the overall effect of small scales is actually to degrade 
the solution for both fields. 

For completeness, we have also decreased the temporal frequency of observation and 
observed that the recovering of the true state was possible provided that the frequency of 
observation, f y , was greater than 3 T^ 1 . 



5 Discussion 



We have derived a quasi-geostrophic model of core dynamics, which aims at describing 
core processes on geomagnetic secular variation timescales. Under the quasi-geostrophic 
assumption, the magnetohydrodynamics takes place in the equatorial plane and is written 
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outside the tangent cylinder. The flow is defined by its zonal velocity and non-zonal 
axial vorticity. The magnetic induction appears through z-averaged quadratic magnetic 
quantities, while we assume that the magnitude of the magnetic field at the core surface is 
smaller than in the core interior. In addition, the equatorial flow is projected on the core- 
mantle boundary. It interacts with the magnetic field at the core surface, through the radial 
component of the magnetic induction equation, in the frozen-flux approximation. That part 
of the model connects the dynamics and the observed secular variation, with the radial 
component of the magnetic field acting as a passive tracer. We have resorted to variational 
data assimilation to construct formally the relationship between model predictions and 
observations. The use of an imperfect observation operator mimics our truncated vision 
of the reality. We have extensively tested the variational data assimilation algorithm 
with twin experiments owing to the non-prohibitive numerical costs of the computations. 
Assimilation was controlled by the initial state and possibly some static model parameters. 

Let us stress some important results from our numerical simulations. In our time- 
dependent framework, we have found that increasing the time window width T always 
improves the solution. In the steady core flow experiment case, that property has proven 
useful in constraining intermediate flow length scales otherwise unconstrained by a sparse 
distribution of observations. Here, the benefit originates from the dynamical relationship 
which exists between successive observations through the radial component of the magnetic 
induction equation. In the torsional oscillation case (an illustration of the dynamical model 
presented in this study), the same property allowed us to retrieve the ^-averaged quadratic 
product of B s (which is only remotely linked to the observed quantity) over the entire 
domain provided that T was large enough for the wave to propagate over (and effectively 
sample) the whole domain. Interestingly, it has not been necessary to include a dissipation 
term in the forward model. We have investigated the sensitivity of the solution to the 
frequency of observation in the presence of observation errors. Adding an extra smoothing 
term to the cost function proved an efficient way to produce analyses with a moderate level 
of complexity. 

From the geophysical point of view, with a typical estimate of the magnetic field 
strength in the core interior, 2 mT, one Alfven time, Ta, amounts to 6 years. We have 
varied the frequency of observation from 2.5 to 20 TJ[ , 3 T^ 1 appearing as a minimum to 
recover the fields, which represents a two-year interval between observations. Therefore, 
we expect that we may be able to resolve properly torsional waves using the last 10 years 
of satellite measur ements, since the most rec ent magnetic field models have a resolution of 
a f raction of year flOlsen and Mandeal. l2008h . 

Pais and Jaultl (l2008h ~ and lGillet et al.l (120091 ) have recently used magnetic field models 



obtained from satellite data in kinematic inversions of quasi-ge ostrophic co r e flow s. Their 
calculated core flows are dominated by a giant retrograde gyre. iGillet et al.l (120091 ) suspect 
that the weaker momentum of the gyre for the period 1960-1980, compared to the period 
1990-2 008, is an artifa c t prod uced by the lesser data quality before 1980. Resorting to a toy 
model, iFournier et al.l (120071 ) have demonstrated that the benefit due to the existence of a 
denser network at the end of an assimilation window, is manifest over a substantial part 
of the window, thanks to the variational data assimilation approach. In other words, the 
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recent high quality of observatory and satellite measurements can be in principle backward 
propagated in order to reassimilate historical data series. It is then possible to imagine 
that the refined series could contribute to a more precise description of both small and 
large scales of the fluid circulation in the core. 

The interplay between magneti c and rotation force s in a two-dimensional model has 
been investigated in other contexts. iTobias et al.l ( 120071 ) have recently studied a local two- 
dimensional /3-plane numerical model to show the impact of a weak large-scale magnetic 
field on the dynamics of the solar tachocline. Instead of using quadratic products of the 
magnetic field as variables, they have written the magnetic field as a function of a unique 
scalar potential A: 



B{s,ip) = V x [A{s,(p)e z 



(4£ 



Then, the magnetic term in the vorticity equation becomes [V x (Ae z )] • V [VfjA] (com- 
pare with the right-hand side term of equation (fl5j) ). and the induction equation is 



-u-VA + S^VlA. 



(49) 



The ansatz (145]) is restrictive, as axial invariance of the magnetic field is assumed. It enables 
the inclusion of magnetic diffusion, the effect of which cannot be rigorously introduced in 
the set of equations (TTTT) to (1T51 . The model given by equations (1451) and (1451 is also 
attractive, since it is still able to describe a variety of physical situations. As an example, 
Diamond et al.l (120051 ) mention the transition from two-dimensional magnetohydrodynamic 
turbulence at small length scales, to turbulence controlled by Rossby wave interactions 
at larger length scales. Solutions of (145|) and (1451 (without the diffusion term) are also 
solutions permitted by our equations (fT5l to (1T51 . Investigating that simplified set of 
equations thus appears as an appealing intermediate step before the actual implementation 
of the less restrictive equations based upon quadratic magnetic quantities. 
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A A first-order variant of the induction equation 



To write equations ( ITT)) to ( TT9)) . we have retained only the zeroth-order part of the 
flow. We may wish to take into account, in these equations, the z-component of the flow 
that enters in the Coriolis term in equation (TT5|) and in the induction equation at the 
core- mantle boundary (equation 01])). 

Thus, in order to ensure incompressibility, we define 

U f( S , V ) = 7(s)[Vx$( s , v )eJ; (50) 

the continuity equation for u E TZ + u\e z yields 7(5) = /i _1 (s). 

The non-zonal vorticity field C, is then defined by £ = V x u^ z , and its vertical 
component is 

( z (s,v) = -V 2 E [h- 1 (s)y(s, V )]. (51) 

Finally, the set of equations (ITT)) to (IT9l becomes 

d t (B 2 S ) = - [u ■ V] (B 2 S ) - 2s- 1 (B 2 S ) u s - 2s- 1 (B 2 ) d v u^ + 2s" 1 (B.BJ d v u a (52) 
d t (Bl) = -[u-V](Bl)-2(B 2 v )d s u s + 2s(B s B ip )d s {s- 1 u ip ), (53) 
d t {B s B v ) = -{u-V](B s B ip ) + s(B 2 )d s (s- 1 u ip )+s- 1 (Bl)d v u s 

- [We ■ u] (B a B v ) . ' (54) 

B Numerical model 

Fields are discretized in radius on a (possibly irregular) staggered grid (see Figure ITU)) . 
s = i s A s (s);i s G [0, N 8 ]; u g {i s + 1/2, j) and r(i s ,j). 1 J r , {B^} and r are calculated on the 
same spatial grid (note that {B^} is not defined on the endpoints). The latitudinal part 
of B r is discretized on a meridian and every grid point is mapped on the CMB from the 
grid point on E, except at the equator (see Figure ITU)) ; thus 9 = ieAg^O); ie G [0, 2N S — 1]. 

The azimuthal part of non-zonal fields, $ and B r is expanded in Fourier series: 

ty(s,<p,t) = 22 [ a m(s,t) cos (rrup) + b m (s , t) sin (rrup)} , (55) 

m=0 

B r (9,(p,t) = ^2 [c m (9,t) cos (mLp) + d m (9,t) sin (rrup)] , (56) 

771=0 

in which the number of azimuthal Fourier mode, m max , is related to the number of equidis- 
tant grid points in longitude, N v : m max = (N v — 1) /2. 

The time step being A t , time is discretized using finite differences t = jA t ; j G [0, N t }. 
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Figure 10: Sketch of a regular staggered radial grid in the equatorial plane, and its projec- 
tion on the CMB. 

Spatial derivatives are computed with a finite difference scheme, except for the longi- 
tudinal derivatives for which we use the Fast Fourier Transform in order to compute them 
in spectral space. 

For the simulations, the cylindrical radius is discretized in N s = 100 grid points (in- 
cluding the boundaries), the CMB in N e = 200 grid points in latitude and N v = 33 (unless 
otherwise specified) grid points in longitude. 

C Steady non-zonal flow model 

We define the streamfunction at the top of the core, \l/ , as = M^\l/, where subscript 
o refers to the outer boundary. Let be the operator which projects a vector from the 
equatorial plane to the top of the core, and let be its transpose. The forward model is 
the radial component of the induction equation at the top of the core, 



B r = -V H -{u NZ B r ), 

= (sin # cos [deVodyBr 



d^ o d B r ] + (cos#) 2 V d v B r . 



(57) 
(58) 



The tangent linear equation is 



5B r = (sin 9 cos Oy 1 [d e V dy6B r 
+ (sin # cos tf)- 1 [d g 5V d v B r 
d t SB r = 5B r . 



d v * d e 6B r ] + (cos ey 2 ^ d^5 B r 

d v 6* dgB r ] + (COS 9)- 2 SVod^Br, 



(59) 



(60) 
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We introduce the adjoint variables iff , and B r for B r and B r respectively. The 
adjoint model is 



B 



B r 



j=0 



(sin 6 cos 6) 1 d ip B rj B r . 



-dl 



(sin 6 cos 9) 1 dgB r . B r . T 



+ (cos6)~ z d v B rj B rj * [> cos 2 9, 



T 



+d'; 1 



fsintfcostfr 1 ^^/ 



(cosfl)" 2 m o B r T 



-dl 



(sin6>cos6>) 1 d ip iff B r T 



(61) 



(62) 



-d t Bj, 



where and dj are the adjoints of d v and dj, and —d t indicates that the integration is 
performed backward in time. 



D Alfven waves model 

The forward equations are: 



d t u g = [s 3 h}~ L d s r, (63) 
d t r = s 3 h{B 2 }d s u g} (64) 
d t B r = -M e s Mlt^u; g d v B ri (65) 

where M total , transforms a one- dimensional zonal vector into a two-dimensional one, by 

zonal " -J 

duplicating it in each meridional plane. Let M^ 1 be its transpose. The forward model is 
completed by the boundary conditions (j4T) and (j4"2"j) . 

We define the adjoint variables wj, r T , {B 2 } T , Bj for u g , r, {-B 2 } , B r respectively. The 
adjoint model is 

d t u T g = dJ[s 3 h{B 2 s }T T }-MZ::;M s e [d^B r Bj} } (66) 
d t r T = dj^h)- 1 ^], (67) 
F T (s) = Y. s * hT J {B 2 s }(d s Lo g ) j+ a c W{B 2 s } } (68) 

3 

d t Bj = -d T v [M^M^Bj] , (69) 

where dj is the adjoint of the operator d s and the term in «p corresponds to the extra 
penalty term in the misfit function (see also equation l2~B"j) ; W = djd s in the experiments. 
In order to enforce its positivity during the optimization phase, {B 2 } is rather written 
{B 2 } = exp [.F(s)], with Fel and F T computed as indicated in equation (168]) . 
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